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We investigate the charge transport property of superconductor (S) /normal metal (N) / fer- 
romagnet insulator (FI) /(normal metal) N’ and S/N/FI/N’/S Josephson junctions on a three- 
dimensional topological insulator surface. We find the asymmetric local density of states (LDOSs) 
in a S/N/FI/N’ junction and show that the N interlayer gives rise to subgap resonant spikes in the 
differential conductance and LDOSs. In a S/N/FI/N’/S junction, the Josephson current shows a 
non-sinusoidal current-phase relation and the N (or N’) interlayer decreases the magnitude of the 
critical current monotonically. 

PACS numbers: 74.45.-|-c,71.10.Pm,74.90.+n 


I. INTRODUCTION 


Three dimensional (3D) topological insulator (TI) is a 
phase of matter with topologically protected Dirac-type 
surface states on their time reversal invariant point [ll- 
[l^ . With coupling to a ferromagnet (F), Dirac fermions 
show many exotic properties such as magnetoelectric ef¬ 
fect |13l4l8l| . By the proximity effect to a superconductor 
(S), the 3D TI surface states may become a topological 
superconductor [l^ . When F and S coexist on 3D TI sur¬ 
faces, it is found that chiral Majorana edge states can be 
generated at the boundary between them [l9l - l2lj| . which 
leads to the formation of zero-biased conductance peak 
(ZBCP) [1^ as experimental signatures!!^ In¬ 

trinsic topological superconductivity has also been found 
in doped 3D TIs, e.g., CuxBi 2 Se 3 

On the other hand, a variety of interesting phe¬ 
nomena about Josephson effect in TI materials have 
been discovered 4ll-l4^. Recently, a non-sinusoidal 
current-phase relation has been reported in the 3D TI 
HgTe junction In the 3D TI heterojunctions 

like Nb/Bii.sSbo.sTei.ySei.s/Nb, the temperature depen¬ 
dence of the critical current is almost linear in most of 
the range [d^. Also, the novel Josephson effect involving 
Majorana fermions has been predicted theoretically (23I - 
[ 2 ^, however, there has been no experimental report yet. 
The rapid development in experiments requires for a the¬ 
oretical approach which can deal with realistic structures 
for Josephson junctions on 3D TI surface. 


In this article, we address how to compose Green’s 
function by wave functions on superconducting 3D TI 
surface. Using the resulting formalism, one can ana¬ 
lyze the spacial dependence of physical quantities, such 
as local density of states (LDOSs) and pair potentials. 
Also, this approach provides an efficient way to cal¬ 
culate Josephson current for realistic junctions on 3D 
TI surfaces. In this work, we consider the S/normal 


metal (N)/ferromagnetic insulator (FI)/N’ junction and 
S/N/FI/N’/S Josephson junction as examples. Since 
making direct contact between F and S regions is not 
easily accessible in actual experiments, the presence of 
N interlayer between S and F is a more realistic setup to 
study Majorana fermions. In the S/N/FI/N’junction, we 
find that the conductance spectra and LDOSs have spikes 
as a function of bias voltage and quasiparticle energy E, 
respectively. The resulting LDOSs shows an asymmetric 
energy dependence around E = 0 . For the S/N/FI/N’/S 
junction, we find that the distance of N interlayer (or N’ 
interlayer) decreases the critical current monotonically. 
The junctions with or without FI show a non-sinusoidal 
current-phase relation at low temperatures. 

The paper is organized as follows: In section II, we in¬ 
troduce our model and construct the Green’s function. In 
section III, we show numerical results for S/N/FI/N’ and 
S/N/FI/N’/S junctions and discuss them. A conclusion 
remark is given in Section IV. 


II. MODEL 


We consider the ballistic S/N/FI/N’ and S/N/FI/N’/S 
junctions which are shown in Fig[T] The system can be 
described by the BdG Hamiltonian [ 2 ^, 


H = 


h{k,c, ky) -P M 


iCTyA 


—iuyA* 


-h*{-k^,-ky) - M* 


( 1 ) 


in (a>^, basis, where h{k,^, ky) = Vf^kyA^ - 

kx^y) /^(G ^ X P J- 0 (x f^) 

the S/N/FI/N’ (S/N/FI/N’/S) junction. ai=x,y,z are 
the Pauli matrices in the spin space and pL is the 
chemical potential. Throughout the paper, we set 
h = 1. The exchange field in F region is M = 
Y.i=x,y,z - ^n(nl)) © (Ln(nl) + - x) for the 
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3D TI Surface 


FIG. 1. Schematics of the system: (a) S/N/FI/N’ and (b) 
S/N/FI/N’/S formed on the surface of a 3D topological insu¬ 
lator. The local density of states can be detected by the STM 
tip. The differential conductance and the suppercurrent can 
be obtained from the leads on the two sides. 


S/N/FI/N’ (S/N/FI/N’/S) junction. The pair poten¬ 
tial A is given by Ao0(—x) for the S/N/FI/N’ junction 
and Ao[0(-x) -I- e~'^Q{x — - Lf — Ln2)] for the 

S/N/FI/N’/S junction, where (j) is the macroscopic su¬ 
perconducting phase. 

In this article, we use a standard formula of tunneling 
spectroscopy 22, 5 ^ as shown in Ref.ll^ to obtain dif¬ 
ferential conductance spectra of the S/N/FI/N’ junction. 
Here, we would like to present the way of constructing 
the retarded Green’s function which has recently been 


applied to relativistic system like Graphene | 54 l l 55 | , and 
ID helical states on TI [^. In our system, the transla¬ 
tional invariance along the y-axis is preserved, thus the 
retarded Green’s function with respect to EqlT] has the 
form G{x,x',y,y') = (x, The re¬ 

tarded Green’s function can be written as IdSMsiL ISSll 


(x, x') = ai-ipi (x)V>f’(x') -I- a2i’i (xj^ij(x') 
-|-a3V’2(x)V’|’(x') -I- a 4 'ip 2 {x)^I{x')> 


with wave vector —ky and h{kx,ky) is given by 
h{kx,ky) = Vf{-kyax - k^ay) - y[Q{-x + dni) + 0(x - 
dni — df)]. For example, in the left S side, the wave 
functions are 

i/i(x) = -b -b ( 5 a) 

Mx) = + 02^6-*'=+" + 52 ^ 6 *'=-", ( 5 b) 

Mx) = C 3 i 3 e-*'=+" + ^3^46*'=-", ( 5 c) 

Mx) = C4i4e*'=-^ -b ^4^36“*'=+“, ( 5 d) 

and 

^i{x') = ( 6 a) 

Mx') = 526 -*'=-"' + d 24 e-*'=+"' + 52546 *'=-"', ( 6 b) 

^3(^0 = 53-836-*'=+"' + ^3546*'=-"', ( 6 c) 

^i{x') = 64.846*'=-"' + ^4-836-*'=+"'. (6d) 


The corresponding wave vectors are represented by k± = 
^J{y± - kl = qe(^h) cos 6»± and qet^h) = 

{y ± \/E'^ — Aq)/?;/. The spinors are given as 

Ai{B'i)=\iu, ±6=*=*®+M, ( 7 a) 

^2(-84)=[i6=*=*®-x, ±u, ie=’=*®- 4 t]^, ( 7 b) 

A^{Bi)=[ie^^^+ u, TM, ±w, *6=’=*®+'u]^, ( 7 c) 

A 4 {B 2 )=\iv, ±e=’=*®-u, =Fe=*=*^-u, iuY', ( 7 d) 

where u and v are given by u(y) = 

Y''( 5 ± - Al)/ 2 E. Other wave functions can 

be found in the Appendix. The coefficients at, bi, 
di and bi can be solved from the boundary con¬ 

dition for relativistic systems. For example, in 
S/N/FI/N’ junction, the boundary conditions are: 
1 pi{x = 0+) = 1 pi{x = 0_), 1 pi{x = dn+) = -ipiXx = dn-), 
ipiix = dn + df+) = X’Xx = dn + df-), and similar to 
other processes. ai=ir ^4 and l3i=i^4 can be determined 
by the boundary conditions of Green’s function 

G'=''(x -b 0 ,x) — G'=“(x — 0,x) = vj ^{ ifza - y ), (8) 


for X > x' and 

G'=s'(x,x') = (diipsix)^^(x') + P2ip4ix)f'[{x') , , 

+I33'ip3{x)lp2 (x') + /34'lp4(x)-lp2 (x')^ 

for X < x'. '(/'i=i~4 (x) are wave functions of Eq. © with 
wave vector ky. 'ipi{2){x) is the wave function for an in¬ 
cident electron-like (hole-like) particle from the left side. 
1/3(4) (2;) is the wave function for the incident electron¬ 
like (hole-like) particles from the right side. '4’i=i~4{x') 
are the wave functions corresponding to the conjugate 
processes under the Hamiltonian 


H = 


h{kx,ky)+M* iayA* 

—iUyA —h*{—kx,—ky) — M 


( 4 ) 


where Ti=x,y,z are the Pauli matrices in the electron-hole 
space. In real materials, the magnitude of the supercon¬ 
ducting gap is much smaller than the chemical potential 
Aq so we can use the quasiclassical approxima¬ 
tion as qe ~ qh and 64. 9 - = 9 . Then one can easily 

obtain the values of ai=ir.,4 and /3i=i.^4, 

01(4) = [ 2 iz;/cos 6 i(M^ - x^)(fi3(i4 - C3C4)]-^C4(3), ( 9 a) 

0 : 2 ( 3 ) = [2i4;/cos6i(M^ - x^)(c3C4 - fi3(i4)]-^(i3(4),(9b) 
/3 i(4 ) = [ 2 w/COs 6 <(m^ - V^){dzd4 - C3C4)]-^C4(3), ( 9 c) 
/32(3) = [2w/cos6*(m^ - X^)(C3C4 - d3d4)]~^d3(4).{9d) 

From the Green’s function, we can obtain the local den¬ 
sity of states for electrons: pe{x,E) and that for holes: 
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Ph{x,E), 

Pe{h){x^E') i?) + p ^ E') ^ (10) 

where the spin-resolved LDOSs are given by 

Pe,t(i) = -\'^'^MG\lf^ 2 )i^,x,E)], (11) 

ky 

PhMl)i^’^) = -\ ^Im[G3%44)(x,x,i?)]. (12) 

ky 

The dc Josephson current is determined by electric 
charge conservation rule 


dtP -l- dxJx -b S' — 0, 


(13) 


where P = -I- Jx = 

S = 2Im[A*4'^\l'-|- — are electric charge den¬ 

sity, electric current and source term, respectively. Af¬ 
ter straightforward calculations following Ref. 4^, M | j we 
find that the total Josephson current is 


Jx = eksT 


ky , UJr ^ 


A sgn(a;„) 

2 v^C)lT^ 


[ai{iuJn) - a2{iuJn )], 


(14) 

where ujn is the Matsubara frequency a;„ = ■nkBT{2n + 
l),(n = 0, ±1,±2..^ Eg. dT^ shows that Furusaki- 
Tsukada’s formula can also be applicable to the 
ballistic Dirac-like electron systems on 3D TI surfaces 
^,5^. It enables us to directly calculate the dc Joseph¬ 


son current in even more complicated or long Josephson 
junctions on 3D TI surface without starting from the en¬ 
ergy levels of Andreev bound states [23|, l59| . 


III. NUMERICAL RESULTS 
A. S/N/FI/N’ junction 



eV/Ao eV/Ao 


FIG. 2. Normalized tunneling conductance as a function of 
bias voltage (eV/Ao) for S/N/FI/N’junctions, (a), (b): = 

0, (c), (d): Ln = i and (e), (f): = 3^. Black curve: 

rriz/p = 1 and red curve: rrix/p = 1. p=^, a/=l, Aq = 0.001 
and Lf = 0.001^ are chosen for all the panels. 


Pe(/i) (ic, E) to that of the electron density of states of 
the bulk normal metal pn at Fermi energy. Here, we 
choose the position in the middle of FI xq = Lf/2 + Ln 
and show the density of states in FiglSl When ^ 0, 
we obtain the subgap peaks again as shown in Figl^Dc ^ 
f). The formation of such peaks can be explained as fol¬ 
lows. We know that the wave vector for electron (hole) is 

k^ = ^{p ± E^jv'j — ky. The condition of forming the 
Andreev bound states in the N layer can be estimated 
from the Bohr-Sommerfeld quantization condition as 


First, we show the conductance CTs (see Appendix) of 
S/N/FI/N’ junction in Figl2j We normalized as by cr„ 
which is the conductance when S is in normal state. We 
only consider the exchange field along z— and a:—axis 
since the magnetization along y—axis does not change 
the conductance [i^. The length of the N layer be¬ 
tween S and FI is denoted by L„. The direct contact 
between S and FI means = 0. For sufficient large 
ixizimx), the normalized conductance has a ZBCP sim¬ 
ilar to that in chiral p-wave superconductor when 
magnetic field is along z-axis as shown in FiglJJa). Also 
we can see from Figl2l(b), ZBCP appears when the mag¬ 
netization is along x—axis. As increases, the sub-gap 
resonant peaks show up (Figsl2](c)^(f)). The number of 
such peaks grows with 

This oscillatory phenomenon can also be seen in 
the local density of states Pejh) (a;, A). We normalize 


edfej-fenlin = ( 15 ) 

which shows that the number of peaks is proportional to 
L„. Similar formation of Andreev bound states was also 
revealed in junctions with ID helical edge states (60j|. 

We also find the asymmetric E dependence of LDOSs 
near the S/FI interface, e.g., pe{xQ,E){ph{xQ,E)) in 
Fig El The asymmetry becomes prominent when mag¬ 
netization is along z axis (FigsEJa), (c) and (e)). We 
know that pe {x, E) and ph [x, E) are symmetric func¬ 
tions of E for chiral p-w ave superconductor when Aq is 
much smaller than p [6ll |. In that case, the time-reversal 
symmetry is already broken in the bulk states of p-wave 
superconductor. On the other hand, the superconductor 
on TI is time-reversal invariant and can not support chi¬ 
ral edge mode without attaching ferromagnet. Therefore, 
we can imagine that the chiral edge mode studied here 
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E/Ao E/Ao 


FIG. 3. Local density of states in the middle of FI in the 
S/N/FI/N’ junctions as a function of energy (E/Aq)'. (a), 
(b): Ln = 0, (c), (d): ^ and (e),(f):L„ = 3^. Solid line 

for electron density of states and dashed line for hole density 
of states. Other parameters are chosen as the same as in 

FigEl 


has a nature similar to Shiba-type bound states 
by magnetic impurity scattering. 




In usual case, where the spin degree of freedom is de¬ 
generate, the emerging Shiba-states still follow the re¬ 
lation Pe(^)(a:, £1) = Pe(h)ix, —E), although the decom¬ 
posed LDOS in each spin sector does not 

satisfy Pe,aix,E) = Pe,a{x,-E). Since Pe{h),crix, E) = 
Pe(h),-a(x,—E) is Satisfied, after summing up each spin 
component, Pe(h){x,E) = p^^^h)^^{x,E) + Pe(h),i{x,E) = 
Pe(h),i{.x,-E) E p^^^h)^-^{x,-E) = p^(^h){x,-E) is satis¬ 
fied. Then, the resulting LDOS is symmetric around 
E = 0. On the other hand, if the spin degeneracy is lifted 
in the superconductor, it is possible that the LDOS be¬ 
comes asymmetric. In the present case, there is a strong 
spin-momentum locking in the superconducting region 
by spin-orbit coupling. Then the asymmetric energy de¬ 
pendence of Pf.(h){x,E) appears near the S/FI interface. 
In recent experiment of scanning tunneling spectroscopy 
(STS), similar asymmetric behavior of LDOSs has been 
observed in ID S/F system]^. We can regard our find¬ 
ing in Fig. [3] as another example of asymmetric LDOSs 
in planar S/F junction which can be detected in STS. 


To see the spacial dependence of the Majorana states 
in such junctions, we show the zero energy density of 
states pe{x,E = 0) throughout the junction. Because 
Pe{x,E = 0) is 0 in both isolated S and FI region, sig¬ 


nificant enhancement of pe{x,E = 0) in S/FI interface 
of S/FI/N junction can be regarded as the experimental 
signature of chiral Majorana fermion. In the S region, 
we can estimate that the characteristic length expressing 
the spatial change oi pe{x, E = 0) is the order of macro¬ 
scopic length scale: This means a sufficient possibility 
to detect the presence of Majorana fermion experimen¬ 
tally by STS, since the manipulation of tip of STS just on 
the the S/N or S/F boundary with high resolution is not 
easy. Also, as seen in FigHJb), even if there is a normal 
layer between S and FI, the enhancement of Pe{x, E = 0) 
in both F and S is not affected. In the N layer between 
S and FI, pe{x,E = 0) is almost constant. In the right 
N layer, we find oscillations oi pe{x, E = 0) on the scale 
of the inverse Fermi momenta. However, this oscillatory 
behavior may be difficult to be detected in actual exper¬ 
iment. 




FIG. 4. Spacial dependence of zero energy states in (a) 
S/FI/N junction and (b) S/N/FI/N’ junction. The width 
of F layer is Lf = 0.001^ and that of the N layer in (b) is 
Ln = C- Other parameters are the same as in Fig[2] The 
scale of the horizontal axis is different in each region. 


B. Josephson effect 

Before discussing the S/N/FI/N’/S junction, let us 
first look at S/N/S junction. Using Ea. dHl) . we plot 
the dc Josephson current in FiglS] It is normalized to 
cEnJ/Aq where Rn is the interface resistance per unit 
area in the normal state. In panel (a), we can see that 
the current-phase relation is non-sinusoidal for short- 
junction in low temperature. This characteristic remains 
in the long-junction, as shown in panel (c). We notice 
that in recent experiment of Nb/3D-HgTe/Nb Josephson 
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FIG. 5. S/N/S Josephson junction: (a) Current-phase rela¬ 
tion and (b) critical current for Ln = 0.01^. (c) and (d) are 
those for Length dependence of Z/„ for T = O.lTc: 

(e) current-phase relation and (f) critical current. Other pa¬ 
rameters are chosen as the same as in Fig[2] 


junctions, the current-phase relation is found to be non- 
sinusoidal 4]^. The experimental condition corresponds 
to low temperature and long-junction in our calculation. 
We find a similar result in that limit as shown in (c). The 
temperature dependence of critical current Jc for short 
and long circumstances are given in panel (b) and (d), 
respectively. We observe that for high temperature, Jc 
is a concave function of T at small while it becomes 
a convex function with large It is also interesting 
to notice that in the large area of low temperature, Jc 
is nearly a linear function of T in both short- and long- 
junction. This result is in good agreement with the recent 
experiments in long Nb/Bii.sSbo.sTei.rSei.a/Nb Joseph¬ 
son junction [46|. In Figsinije) and (f), we plot the length 
dependence of Josephson current. 


We now consider S/N/FI/N’/S Josephson junctions. 
The length of N layer on the two side of FI is denoted 
as Lni and L„ 2 - When Lni and Ln 2 is on the super¬ 
conducting coherence length scale, the junctions become 
long-junctions. The influence of N layer between FI and S 
is shown in Figl6l From FigslU^a) and (b), we can see that 
the current-phase relation still retains the non-sinusoidal 
shape for different values of L„i and L „2 in low temper¬ 
ature limit. Throughout our study, we have not found 
the sawtooth behavior of current-phase relation involving 


magnetization in the long junction and low temperature 
limit, as shown in FiglSjc). This is because the mag¬ 
netization makes the Andreev bound states gapped for 


most values of 


The derivative of energy disper¬ 


sion which creates Josephson current will be a smoother 
function of phase than that in S/N/S junction. For the 
temperature dependence of the critical current, we can 
see that it behaves qualitatively different in low temper¬ 
ature limit for and rux as shown in FigslU^c) and (d), 
respectively. For rriz case, the critical current Jc satu¬ 
rates at a constant value, which has been revealed by the 
previous work 2g. However, for the nix case^h shows 


a Kulik-Ome’lyanchuk type of critical current [62| which 
has linear low-temperature behavior. We interpret it as 
a result of the enhanced zero-energy LDOSs for nix mag¬ 
netization as illustrated in Figsl31[b)(d)(f). In the high 
temperature limit, it is shown that for both and nix 
cases, Jc is a concave function while it crosses over to 
a convex function with increasing Lni (or Ln 2 )- This 
behavior is similar to the S/N/S junction. FigureslU^e) 
and (f) represent the critical current as a function of the 
length Lni and Ln 2 , for different direction of magneti¬ 
zation. It is worth noting that, although the interlayer 



FIG. 6. S/N/FI/N’/S junction: current-phase relation for 3 
cases of the N layer length Lni and Ln 2 for (a) niz/fi = 1 
and (b) ~ 1. (c)(d) Temperature dependence of criti¬ 

cal current corresponding to (a) and (b), respectively. (e)(f) 
Gritical Josephson current as a function of Lni and Ln 2 ■ The 
temperature is chosen as T — O.lTc. Other parameters are 
chosen as the same as in Fig[5] 
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N in the S/N/FI/N’ junction could generate resonant 
spikes in the transport phenomena, e.g., spikes in Figsl^ 
and [21 we find no oscillatory behavior in either current- 
phase relation or critical current as a function of length 
N (or N’). The critical current decreases monotonically 
with the length -I- L„ 2 - 


with 

knM2) = ( 20 ) 

^n.3(4) = (a^ - EY/v'j -kl = ±K- (21) 

For the FI interlayer, we find that 


IV. SUMMARY 

In summary, we theoretically studied the S/N/FI/N’, 
S/N/S and S/N/FI/N’/S junctions on the surface of 3D 
topological insulator. We have constructed a formula to 
obtain Green’s function. The conductance spectra and 
local density of states in S/N/FI/N’ junction show reso¬ 
nant spikes due to the Andreev bound states. The cal¬ 
culated current phase relation and temperature depen¬ 
dence of critical current in the Josephson junctions are 
consistent with recent experiments in S/N/S junction. 
We have also calculated current phase relation and tem¬ 
perature dependence of critical current in S/N/FI/N’/S 
junction. The non-sinusoidal current phase relation can 
be expected for short junctions. We hope the obtained 
results will be confirmed by experiments in the near fu¬ 
ture. 


A=1 


A=1 


( 22 ) 

(23) 


where 
Fi = 

F 2 = 

h = 


Fa = 


ivfk{ {vfky + mx) ,E — , 

f 

E + mz,—ivfk^ + {vfky + rrix) ,0,0 , 

0,0,—ivfk^ + {vfky — mx),E + mz , 
f 

0,0,E - mz,ivfk^ + {vfky - ruj;) , 


(24a) 

(24b) 

(24c) 

(24d) 


with 
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APPENDIX: WAVE FUNCTIONS 

The wave functions in the N interlayer are 

4 

= (16) 

A^l 
4 

= (17) 

A=1 

where 


k{ = -^i\JE'^ - ml - {vfky + 

( 25 a) 

ki = ‘^i\jF^ - ml - [vfky + mxf', 

( 25 b) 

H = ^2\/ E"^ - ml - {vfky - m^f, 

( 25 c) 

ki = -^ 2 \lE'^ - ml - {vfky - mxf. 

( 25 d) 

and <ri(2) = sgn (vfky ± mx)- The wave functions 
N’ region of S/N/FI/N’ are 

in the 

V'i(x) = 

( 26 a) 

Mx) = C2C2e-*'=" " + d2Gie*'="", 

( 26 b) 

Mx) = 46 -*'=-“ + 03(726-*'=" " + 

( 26 c) 

Mx) = + 64(726-*'="’’=. 

( 26 d) 

and 


Mx) = 5 iDi 6 *'=""' + Jii 72 e-*'=""', 

( 27 a) 

Mx) = £2^26-*'=""' + ^2^16*'=""', 

( 27 b) 

Mx') = 1736-*'=""' + d3D26-*'=7"' + 63 Di6*'=""', 

( 27 c) 

Mx') = .D 46 *'= 7 "' + d 4 Di 6 *'=""' + 64 .D 26 -*'= 7 "'. 

( 27 d) 


j, with 

^i{2) = [I’f {'i-kn,i{2) + ky) , fj, + E,0,0] , (18) 

^3(4) = [0,0, Vf {ikn,3{A) - ky) , -+ E]^ , (19) 


kt = cos6»±. 


(28) 
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The spinors are given by 


Ci(i)3)=[*,±e±*"”,0,0]^, 

(29a) 

C2(£>4)=[0,0,±l,ie±*®"]^, 

(29b) 

C3pi)=[*e±<,Tl,0,0f, 

(29c) 

C'4(£>2)=[0,0,Te±*®",^]^. 

(29d) 

Also, the conductance can be given as 


0 's = o-Q J dkyRe 1 + ^ ^ - Ihf 

(30) 


where (Tq is a constant parameter determined by the ge¬ 
ometry of junctions. 
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